1. Field of the Invention
The present invention relates to statistical methods for reconstructing a polyenergetic X-ray computed tomography image and image reconstructor apparatus and, in particular, to methods and reconstructor apparatus which reconstruct such images from a single X-ray CT scan having a polyenergetic source spectrum.
2. Background Art
X-ray computed or computerized tomography (i.e. CT) provides structural information about tissue anatomy. Its strength lies in the fact that it can provide xe2x80x9cslicexe2x80x9d images, taken through a three-dimensional volume with enhanced contrast and reduced structure noise relative to projection radiography.
FIG. 1 illustrates a simple CT system. An X-ray source is collimated and its rays are scanned through the plane of interest. The intensity of the X-ray photons is diminished by tissue attenuation. A detector measures the photon flux that emerges from the object. This procedure is repeated at sufficiently close angular samples over 180xc2x0 or 360xc2x0. The data from different projections are organized with the projection angles on one axis and the projection bins (radial distance) on the other. This array is referred to as the sinogram, because the sinogram of a single point traces a sinusoidal wave. Reconstruction techniques have the goal of estimating the attenuation map of the object that gave rise to the measured sinogram.
FIGS. 2a-2c illustrate the evolution of CT geometries. FIG. 2a is a parallel-beam (single ray) arrangement, much like what was found in a first-generation CT scanner. The major drawback of this arrangement is long scan time, since the source detector arrangement has to be translated and rotated.
The fan-beam geometry of FIG. 2b reduces the scan time to a fraction of a second by eliminating the need for translation. By translating the patient table as the source detector arrangement rotates, one gets an effective helical path around the object leading to increased exposure volume and three-dimensional imaging.
The latest CT geometry is the cone-beam arrangement, shown in FIG. 2c. It further reduces scan time by providing three-dimensional information in one rotation. It is most efficient in its usage of the X-ray tube, but it suffers from high scatter (xe2x89xa740%). It is also the most challenging in terms of reconstruction algorithm implementation.
Two dominant effects, both a function of the X-ray source spectrum, govern tissue attenuation. At the lower energies of interest in the diagnostic region, the photoelectric effect dominates. At higher energies, Compton scattering is the most significant source of tissue attenuation.
The linear attenuation coefficient xcexc(x,y,z,E) characterizes the overall attenuation property of tissue. It depends on the spatial coordinates and the beam energy, and has units of inverse distance. For a ray of infinitesimal width, the mean photon flux detected along a particular projection line Li is given by:
E[Yi]=∫Ii(E)exe2x88x92∫Lixcexc(x,y,z,E)dldExe2x80x83xe2x80x83(1)
where the integral in the exponent is taken over the line Li and Ii(E) incorporates the energy dependence of the incident ray and detector sensitivity. The goal of any CT algorithm is to reconstruct the attenuation map xcexc from the measured data [Y1, . . . , YN] where N is the number of rays measured.
Filtered back projection (FBP) is the standard reconstruction technique for X-ray CT. It is an analytic technique based on the Fourier slice theorem.
Use of the FFT in the filtering step of FBP renders the algorithm quite fast. Moreover, its properties are well understood. However, because it ignores the noise statistics of the data, it results in biased estimators. It also suffers from streak artifacts when imaging objects with metallic implants or other high-density structures.
In reality, the attenuation coefficient xcexc is energy dependent and the X-ray beam is polyenergetic. Lower energy X-rays are preferentially attenuated. FIG. 7 shows the energy dependence of the attenuation coefficients of water (density 1.0 gm/cm3) and bone (at density 1.92 gm/cm3). A hard X-ray beam is one with higher average energy. Beam hardening is a process whereby the average energy of the X-ray beam increases as the beam propagates through a material. This increase in average energy is a direct consequence of the energy dependence of the attenuation coefficient.
With a polyenergetic source, the expected detected photon flux along path Li is given by (1). If one were to ignore the energy dependence of the measurements and simply apply FBP to the log processed data, some attenuation map {circumflex over (xcexc)} would be reconstructed that is indirectly related to the source spectrum and object attenuation properties.
Beam hardening leads to several disturbing artifacts in image reconstruction. FIGS. 8 and 9 show the effect of beam hardening on the line integral in bone and water. In the monoenergetic case, the line integral increases linearly with thickness. With a polyenergetic beam, the soft tissue line integral departs slightly from the linear behavior. The effect is more pronounced for high Z (atomic number) tissue such as bone.
This non-linear behavior generally leads to a reduction in the attenuation coefficient. In bone, beam hardening can cause reductions of up to 10%. Thick bones also generate dark streaks. In soft tissue, the values are depressed in a non-uniform manner, leading to what has been termed the xe2x80x9ccuppingxe2x80x9d effect. In addition, bone areas can xe2x80x9cspill overxe2x80x9d into soft tissue, leading to a perceived increase in the attenuation coefficient.
Because of SNR considerations, monoenergetic X-ray scanning is not practical. Beam hardening correction methods are therefore necessary for reconstructing artifact-free attenuation coefficient images. An ideal reconstruction method would be quantitatively accurate and portable to different scanning geometries. It would somehow reconstruct xcexc(x,y,E), retaining the energy dependence of the attenuation process. This is difficult, if not impossible, to achieve with a single source spectrum. A more realistic goal is to remove or reduce the beam hardening artifacts by compensating for the energy dependence in the data.
There are a wide variety of schemes for beam hardening artifact reduction. Existing methods fall into three categories: dual-energy imaging, preprocessing of projection data and post-processing of the reconstructed image.
Dual-energy imaging has been described as the most theoretically elegant approach to eliminate beam hardening artifacts. The approach is based on expressing the spectral dependence of the attenuation coefficient as a linear combination of two basis function, scaled by constants independent of energy. The two basis functions are intended to model the photo-electric effect and Compton scattering. This technique provides complete energy dependence information for CT imaging. An attenuation coefficient image can, in principle, be presented at any energy, free from beam hardening artifacts. The method""s major drawback is the requirement for two independent energy measurements. This has inhibited its use in clinical applications, despite the potential diagnostic benefit of energy information. Recently, some work has been presented on the use of multi-energy X-ray CT for imaging small animals. For that particular application, the CT scanner was custom built with an energy-selective detector arrangement.
Commercial beam hardening correction methods usually involve both pre-processing and post-processing, and are often implemented with a parallel or fan-beam geometry in mind. They also make the assumption that the object consists of soft tissue (water-like) and bone (high Z). Recently, these methods were generalized to three base materials and cone-beam geometry.
Pre-processing works well when the object consists of homogeneous soft tissue. Artifacts caused by high Z materials such as bone mandate the use of post-processing techniques to produce acceptable images.
The attenuation coefficient of some material k is modeled as the product of the energy-dependent mass attenuation coefficient mk(E) (cm2/g) and the energy-independent density xcfx81(x,y) (g/cm3) of the tissue. Expressed mathematically,                               μ          ⁡                      (                          x              ,              y              ,              E                        )                          =                              ∑                          k              =              1                        K                    ⁢                                                    m                k                            ⁡                              (                E                )                                      ⁢                          xe2x80x83                        ⁢                                          ρ                k                            ⁡                              (                                  x                  ,                  y                                )                                      ⁢                          xe2x80x83                        ⁢                                          r                k                            ⁡                              (                                  x                  ,                  y                                )                                                                        (        2        )            
where k is the number of tissue types in the object and rk(x,y)=1 if (x,y) xcex5 tissue k and rk(x,y)=0 otherwise.
For the classical pre-processing approach, the object is assumed to consist of a single tissue type (K=1) and to have energy dependence similar to that of water, i.e.,                               μ          ⁡                      (                          x              ,              y              ,              E                        )                          =                  xe2x80x83                ⁢                                            m              1                        ⁡                          (              E              )                                ⁢                      xe2x80x83                    ⁢                                    ρ              1                        ⁡                          (                              x                ,                y                            )                                                                    xe2x80x83                ⁢                  (          3          )                                                  ≈                      xe2x80x83                    ⁢                                                    m                w                            ⁡                              (                E                )                                      ⁢                          xe2x80x83                        ⁢                                          ρ                soft                            ⁡                              (                                  x                  ,                  y                                )                                                    ⁢                  xe2x80x83                                              xe2x80x83                ⁢                  (          4          )                    
where mw(E) is the mass attenuation coefficient of water and xcfx81soft is the effective soft tissue density. One can rewrite (1) as follows:                               E          ⁡                      [                          Y              i                        ]                          =                  xe2x80x83                ⁢                  ∫                                                    I                i                            ⁡                              (                E                )                                      ⁢                          xe2x80x83                        ⁢                          ⅇ                              -                                  ∫                                                            L                      i                                        ⁢                                          μ                      ⁡                                              (                                                  x                          ,                          y                          ,                          E                                                )                                                              ⁢                                          ⅆ                      t                                                                                            ⁢                          ⅆ              E                                                                    xe2x80x83                ⁢                  (          5          )                                                  =                      xe2x80x83                    ⁢                                    ∫                              L                i                                      ⁢                                                            I                  i                                ⁡                                  (                  E                  )                                            ⁢                              xe2x80x83                            ⁢                              ⅇ                                                      -                                                                  m                        w                                            ⁡                                              (                        E                        )                                                                              ⁢                                                            ∫                                              L                        i                                                              ⁢                                                                  ρ                        soft                                            ⁡                                              (                                                  x                          ,                          y                                                )                                                                                                        ⁢                              ⅆ                E                                                    ⁢                  xe2x80x83                                    xe2x80x83                                =                  xe2x80x83                ⁢                              ∫                          L              i                                ⁢                                                    I                i                            ⁡                              (                E                )                                      ⁢                          xe2x80x83                        ⁢                          ⅇ                                                -                                                            m                      w                                        ⁡                                          (                      E                      )                                                                      ⁢                                  T                  i                                                      ⁢                          ⅆ              E                                                                    xe2x80x83                ⁢                  (          6          )                                        =                  xe2x80x83                ⁢                              F            i                    ⁡                      (                          T              i                        )                                                        xe2x80x83                ⁢                  (          7          )                    
where Ti is the line integral of the density along path Li. Each function Fi(Ti) i=1, . . . , N is 1xe2x88x921 and monotone decreasing and hence invertible. The goal of the pre-processing method is to estimate {Ti}i=1N and from that reconstruct (using FBP) an estimate {circumflex over (xcfx81)}(x,y) of the energy-independent density xcfx81soft. In other words, {circumflex over (xcfx81)}(x,y)=FBP {{circumflex over (T)}i} where {circumflex over (T)}i=Fixe2x88x921(Yi). This pre-processing approach is inaccurate when bone is present, but is often the first step in a post-processing bone correction algorithm.
Post-processing techniques first pre-process and reconstruct the data for soft tissue correction, as explained above. The resulting effective density image is then segmented into bone and soft tissue. This classification enables one to estimate the contributions of soft tissue and bone to the line integrals. These estimates are used to correct for non-linear effects in the line integrals. The final artifact-free image is produced using FBP and displays density values independent of energy according to the following relationship:
{circumflex over (xcfx81)}(x,y)={circumflex over (xcfx81)}soft(x,y)+xcexo{circumflex over (xcfx81)}bone(x,y)xe2x80x83xe2x80x83(8)
where xcexo is some constant independent of energy that maintains image contrast.
Although post-processing accomplishes its goal of eliminating energy dependence, it suffers from quantitative inaccuracy. As explainer elsewhere, the parameter xcexo is somewhat heuristically estimated. In addition, applying post-processing to 3-D or cone-beam geometries can be computationally expensive. Another beam hardening correction of interest is known. This algorithm is iterative (but not statistical). At each pixel, it assumes that the attenuation coefficient is a linear combination of the known attenuation coefficients of two base materials, and it iteratively determines the volume fractions of the base materials. The algorithm depends on an empirically-determined estimate of effective X-ray spectrum of the scanner. The main limitations of this approach is that the spectrum estimate captures the imaging characteristics for a small FOV only and that prior knowledge of the base materials at each pixel is necessary.
Many of the inherent shortcomings of FBP are adequately (and naturally) compensated for by statistical methods.
Statistical methods are a subclass of iterative techniques, although the two terms are often used interchangeably in the literature. The broader class of iterative reconstruction techniques includes non-statistical methods such as the Algebraic Reconstruction Technique (ART) which casts the problem as an algebraic system of equations. Successive substitution methods, such as Joseph and Spital""s beam-hardening correction algorithm, are also iterative but not statistical. Hence, statistical methods are iterative, but the opposite is not necessarily true.
Statistical techniques have several attractive features. They account for the statistics of the data in the reconstruction process, and therefore lead to more accurate estimates with lower bias and variance. This is especially important in the low-SNR case, where deterministic methods suffer from severe bias. Moreover, statistical methods can be well suited for arbitrary geometries and situations with truncated data. They easily incorporate the system geometry, detector response, object constraints and any prior knowledge. Their main drawback (when compared to FBP) is their longer computational times.
Statistical reconstruction for monoenergetic CT was shown to outperform FBP in metal artifact reduction, in limited-angle tomography and to have lower bias-noise curves.
All single-scan statistical X-ray reconstruction algorithms assume (either implicitly or explicitly) monoenergetic X-ray beams, and thus do not deal with the issue of beam hardening artifacts. The future potential of statistical methods to correct for beam-hardening was speculated by Lange et al. in 1987, but no methods have been proposed for realizing this potential. One exception is the area of dual-energy CT. Dual-energy systems operate based on the principle that the attenuation coefficient can be expressed as a linear combination of two energy basis functions and are capable of providing density images independent of energy.
For X-ray CT images, with typical sizes of 512xc3x97512 pixels or larger, statistical methods require very long computational times. This has hindered their widespread use.
The article of Yan et al. in Jan. 2000 IEEE Trans. Med. Im. uses a polyenergetic source spectrum, but it is not statistical and it does not have any regularization. There is no mathematical evidence that their algorithm will converge.
In general, all previous algorithms for regularized statistical image reconstruction of X-ray CT images from a single X-ray CT scan have been either:
1) based explicitly on a monoenergetic source assumption, or
2) based implicitly on such an assumption in that the polyenergetic spectrum and resulting beam hardening effects were disregarded.
An object of the present invention is to provide a method for reconstructing a polyenergetic X-ray computed tomography image and an image reconstructor apparatus, both of which utilize a statistical algorithm which explicitly accounts for a polyenergetic source spectrum and resulting beam hardening effects.
Another object of the present invention is to provide a method for reconstructing a polyenergetic X-ray computed tomography image and an image reconstructor apparatus, both of which utilize a statistical algorithm which is portable to different scanner geometries.
In carrying out the above objects and other objects of the present invention, a method for statistically reconstructing a polyenergetic X-ray computed tomography image to obtain a corrected image is provided. The method includes providing a computed tomography initial image produced by a single X-ray CT scan having a polyenergetic source spectrum. The initial image has components of materials which cause beam hardening artifacts. The method also includes separating the initial image into different sections to obtain a segmented image and calculating a series of intermediate corrected images based on the segmented image utilizing a statistical algorithm which accounts for the polyenergetic source spectrum and which converges to obtain a final corrected image which has significantly reduced beam hardening artifacts.
The step of calculating may include the steps of calculating a gradient of a cost function having an argument and utilizing the gradient to minimize the cost function with respect to its argument. The step of calculating the gradient may include the step of back projecting. The cost function preferably has a regularizing penalty term.
The step of calculating the gradient may include the step of calculating thicknesses of the components. The step of calculating thicknesses may include the step of reprojecting the segmented image.
The step of calculating the gradient may include the step of calculating means of data along paths and gradients based on the thicknesses of the components.
The argument may be density of the materials at each image voxel.
The method may further include calibrating the spectrum of the X-ray CT scan.
The method may further include displaying the final corrected image.
The step of calculating the gradient may include the step of utilizing ordered subsets to accelerate convergence of the algorithm.
Further in carrying out the above objects and other objects of the present invention, an image reconstructor apparatus for statistically reconstructing a polyenergetic X-ray computed tomography image to obtain a corrected image is provided. The apparatus includes means for providing a computed tomography initial image produced by a single X-ray CT scan having a polyenergetic source spectrum wherein the initial image has components of materials which cause beam hardening artifacts. The apparatus further includes means for separating the initial image into different sections to obtain a segmented image and means for calculating a series of intermediate corrected images based on the segmented image utilizing a statistical algorithm which accounts for the polyenergetic source spectrum and which converges to obtain a final corrected image which has significantly reduced beam hardening artifacts.
The means for calculating may include means for calculating a gradient of a cost function having an argument and means for utilizing the gradient to minimize the cost function with respect to its argument. The cost function preferably has a regularizing penalty term.
The means for calculating the gradient may include means for back projecting.
The means for calculating the gradient may include means for calculating thicknesses of the components.
The means for calculating thicknesses may include means for reprojecting the segmented image.
The means for calculating the gradient may include means for calculating means of data along paths and gradients based on the thicknesses of the components.
The argument may be density of the materials at each image voxel.
The means for calculating the gradient may include means for utilizing ordered subsets to accelerate convergence of the algorithm.
The above objects and other objects, features, and advantages of the present invention are readily apparent from the following detailed description of the best mode for carrying out the invention when taken in connection with the accompanying drawings.